/*************************************************************
** THROBIT2: Threshold observability bivariate probit model
** Monte Carlo Comparison of Throbit with Boolean Probit
** Skewed to have very few successes
** In this case, one of the risks (the second) only has the common covariate
** (c) 2003 by Sanford C. Gordon
**
**
** Department of Politics
** new York University
** Purpose: Determine sample properties of partially observed probit with thresholds.
** output file: throbit1.out
*************************************************************/
rndseed 33333;
closeall; 
cls; 

/* Call to the Maxlik library */
library maxlik;
/* Resets global variables to default values */
maxset;

/* Global Parameter Settings */

_max_algorithm = 4;
	                 	@ Optimization Method           @
                       	@ 1 = Steepest Descent          @
                       	@ 2 = BFGS                      @
                       	@ 3 = DFP                       @
                       	@ 4 = Newton-Raphson            @
                       	@ 5 = BHHH                      @
                       	@ 6 = PRCG                      @

_max_CovPar = 1;
						@ Covariance Matrix Method		@
						@ 0 = Not Computed			    @
						@ 1 = Numerical Hessian			@
						@ 2 = Quasi-maximum likelihood  @

_max_LineSearch = 2;
						@ Line Search Method			@
						@ 1 = unit step length			@
						@ 2 = cubic or quadratic		@
						@ 3 = step halving			    @
						@ 4 = Brent's method			@
						@ 5 = BHHH method				@
_max_maxiters = 100;

/*__output=0;*/
/*Create symbolic infinity*/
infty=2^10000;


/*******************************************************************************
** Throbit Liklihood Procedure                                                **
*******************************************************************************/

proc throbiti(parms,dta);
    local obs,y,y1,y2,x1,x2,beta1,beta2,tau1,tau2,idx1,idx2,p0,py1,py2,pya,llike;
    /* Data */
    obs = rows(dta);
    y   = dta[.,1];
    y1  = dta[.,2];
    y2  = dta[.,3];
    x1  = ones(obs,1)~dta[.,4:5];
    x2  = ones(obs,1)~dta[.,4];
    /* parameters */
    beta1 = parms[1:3];
    beta2 = parms[4:5];
    tau1  = exp(parms[6]);
    tau2  = exp(parms[7]);
    idx1= x1*beta1;
    idx2=x2*beta2;
    p0 = cdfn(-idx1).*cdfn(-idx2); /* This is the probability of a no */
    py1 = cdfn(idx1-tau1).*cdfn(tau2-idx2); /* This is the probability of a yes via mechanism 1 */
    py2 = cdfn(idx2-tau2).*cdfn(tau1-idx1); /* Probability of a yes via mechanism 2 */
    pya = 1-p0-py1-py2; /* Probability of an ambiguous yes */
    /* This is the likelihood function */
    llike = (1-y).*ln(p0)+y.*y1.*ln(py1)+y.*y2.*ln(py2)+y.*(1-y1).*(1-y2).*ln(pya);
    retp(llike);
endp;

/************************************************************************
** Boolean Probit Liklihood Procedure                                  **
************************************************************************/

proc boolprob(parms,dta);
    local y,llike,obs,x1,x2,beta1,beta2,p1,p2,yesprob,noprob;
    obs=rows(dta);
    y=dta[.,1];
    x1  = ones(obs,1)~dta[.,4:5];
    x2  = ones(obs,1)~dta[.,4];
    beta1=parms[1:3];
    beta2=parms[4:5];
    p1 = cdfn(x1*beta1);
    p2 = cdfn(x2*beta2);
    yesprob = 1-(1-p1).*(1-p2);  /* Probability of a yes */
    noprob  = 1-yesprob; /* Probability of a no */
    /* This is the likelihood function */
    llike = y.*ln(yesprob)+(1-y).*ln(noprob);
    retp(llike);
endp;

/* Parameter Values Fixed for all simulations */
beta1=1|1|1;
beta2=1|1;
tau1 = 1;
tau2 = 0.001;
mciter=1000;    /* Number of Monte Carlo simulations */
xoffset = 3.77; /* When drawing Xs, substract this amount to control the outcome balance */
obsformc = 500|1000|5000; /* Various sample sizes for different experiments */

/*****************************************************************************************
** Monte Carlo Routine: Compare throbit and Boolean probit, for samples of varying sizes**
*****************************************************************************************/


j=1; do until j > rows(obsformc);
    N = obsformc[j];
    /* Generate Xs, which are fixed in repeated samples */
    xcommon = rndn(N,1)-3;
    x12 = rndn(N,1)+0.7;
    x1 = ones(N,1)~xcommon~x12;
    x2 = ones(N,1)~xcommon;
    i=1; do until i >mciter;
    cls; 
    print "N=" n ", i=" i;
    /* Generate observed outcomes */
        epsilon1 = rndn(N,1);
        epsilon2 = rndn(N,1);
        zstar1  = x1*beta1+epsilon1;
        zstar2  = x2*beta2+epsilon2;
        y = maxc((zstar1~zstar2)').>0;
        y1= (zstar1.>tau1).*(zstar2.<tau2);
        y2= (zstar2.>tau2).*(zstar1.<tau1);
        let label = y y1 y2 xcommon x12;
    /* Save the data */
        create f0= "c:/gauss50/research/throbit/throbit1" with ^label,0,8;
        call writer(f0,y~y1~y2~xcommon~x12);
        call close(f0);
    
        dta = "c:/gauss50/research/throbit/throbit1"; /* Tells Maxlik where to look for the data */
        vars = getname(dta);
        
    
        /* First, run throbit */
        _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c"|"lntau1"|"lntau2";
        theta0th=0|0|0|0|0|0|0;
        {thetath,fth,gth,covth,retth}=maxprt(maxlik(dta,vars,&throbiti,theta0th));
    
        /* Next, run Boolean probit, using starting values randomly perturbed around zero 
           (Gauss sometimes has trouble with a zero vector) */
        _max_parnames = "beta10"|"beta1c"|"beta12"|"beta20"|"beta2c";
        theta0bp=0.01*rndn(5,1);
        {thetabp,fbp,gbp,covbp,retbp}=maxprt(maxlik(dta,vars,&boolprob,theta0bp));
    
        
        /* Save the results */
        miss5 = {. . . . .};
        miss7 = {. . . . . . .};
        if i == 1;
            meany = meanc(y~y1~y2)';
            /* Results for Boolean probit */
            resultbp = thetabp';
            convrgbp = retbp;
            llikebp = fbp*N;
            truthbp = beta1|beta2;
            if retbp == 0;
                msematbp = (thetabp-truthbp)*(thetabp-truthbp)';
            else;
                msematbp = zeros(5,5);
            endif;
            if ismiss(covbp) == 0;
                serrbp = sqrt(diag(covbp))'; 
            else;
                serrbp = miss5;
            endif;
            /* Results for throbit */
            resultth = thetath';
            convrgth = retth;
            lliketh = fth*N;
            truthth = truthbp|ln(tau1)|ln(tau2);
            /* Calculates Mean Squared Error Matrices; #2 is for when Boolean probit fails */
            if retbp == 0;
                msematth = (thetath-truthth)*(thetath-truthth)';
                msemtth2 = zeros(7,7);
            else;
                msemtth2 = (thetath-truthth)*(thetath-truthth)';
                msematth = zeros(7,7);
            endif;
            if ismiss(covth) == 0;
                serrth = sqrt(diag(covth))'; 
            else;
                serrth = miss7;
            endif;
        else;
            meany = meany|meanc(y~y1~y2)';
            resultbp = resultbp|thetabp';
            convrgbp = convrgbp|retbp;
            llikebp = llikebp|fbp*n;
            truthbp = beta1|beta2;
            if retbp ==0; 
                msematbp = msematbp+(thetabp-truthbp)*(thetabp-truthbp)';         
            endif;             
            if ismiss(covbp)==0;
                serrbp=serrbp|sqrt(diag(covbp))';
            else;
                serrbp = serrbp|miss5;
            endif;
    
            resultth = resultth|thetath';
            convrgth = convrgth|retth;
            lliketh = lliketh|fth*n;
            truthth = truthbp|ln(tau1)|ln(tau2);
            if retbp == 0;
                msematth = msematth+(thetath-truthth)*(thetath-truthth)';
            else;   
                msemtth2 = msemtth2+(thetath-truthth)*(thetath-truthth)';
            endif;
            if ismiss(covth)==0;
                serrth=serrth|sqrt(diag(covth))';
            else;
                serrth = serrth|miss7;
            endif;
        endif;
        i=i+1;
    endo;    

    sebp500 = meanc(packr(serrbp));
    seth500 = meanc(packr(serrth));
    
    keepbp = convrgbp .== 0;  /* This is a vector of zeros and ones; ones correspond to instances where
                                BP had a hard time */
    numkeep = sumc(keepbp); /* This is the number of simulations where BP found an answer */
    numlose = mciter-numkeep; /* This is the number of simulations where BP failed -- "Hard" Cases */

    rmsebp = sqrt(diag(msematbp))/numkeep; /* Calculate root mean squared error for successful cases of BP */
    rmseth = sqrt(diag(msematth))/numkeep; /* Calculate root mean squared error for throbit estimates 
                                          corresponding to cases in which BP succeeded */
    if numlose == 0;
        rmseth2 = zeros(7,1);
    else;
        rmseth2= sqrt(diag(msemtth2))/numlose; /* Calculate rmse for throbit when BP failed */
    endif;
    
    resbpk = delif(resultbp,1-keepbp); /*Results for successful Boolean probit */
    resthk = delif(resultth,1-keepbp); /*Results for throbit when BP was successful */
    resthl = delif(resultth,keepbp);   /*Results for throbit when BP failed */
    
    /* Print output table */ 
    amper4= "&"$|"&"$|"&"$|"&";
    slashes4  = "\\\\"$|"\\\\"$|"\\\\"$|"\\\\";
    caption = "Truth"|"Mean Boolean"|"Mean Throbit"|"Rel. Eff.";
    meanbp = meanc(resbpk);
    meanth = meanc(resthk);
    releff = rmsebp./(rmseth[1:5]);
    b10 = truthbp[1]|meanbp[1]|meanth[1]|releff[1];
    b1c = truthbp[2]|meanbp[2]|meanth[2]|releff[2];
    b12 = truthbp[3]|meanbp[3]|meanth[3]|releff[3];
    b20 = truthbp[4]|meanbp[4]|meanth[4]|releff[4]; 
    b2c = truthbp[5]|meanbp[5]|meanth[5]|releff[5];
    
    lnt1 = meanth[6];
    lnt2 = meanth[7];
    lnt1t = ln(tau1);
    lnt2t = ln(tau2);
    
    b10s = "" $+ftocv(b10,0,2);
    b1cs = "" $+ftocv(b1c,0,2);
    b12s = "" $+ftocv(b12,0,2);
    b20s = "" $+ftocv(b20,0,2);
    b2cs = "" $+ftocv(b2c,0,2);
    
    lnt1s = "" $+ftocv(lnt1t,0,2)|"--"|"" $+ftocv(lnt1,0,2)|"--";
    lnt2s = "" $+ftocv(lnt2t,0,2)|"--"|"" $+ftocv(lnt2,0,2)|"--";
    obshead = obsformc[j];
    
       
    resmat = caption$~amper4$~b10s$~amper4$~b1cs$~amper4$~b12s$~amper4$~b20s$~amper4$~b2cs$~amper4$~lnt1s$~amper4$~lnt2s$~slashes4;
    if j == 1;
        output file = "c:/gauss50/research/throbit/mc3.out" reset;
        print"\\begin{table}";
        print"\\caption{\\label{mc3} Monte Carlo Comparison of Boolean Probit and Throbit, Minimal Identifying Conditions}";
        print"\\begin{center}";
        print"\\begin{tabular}{lccccccc}";
        print"\\hline \\hline";
        print"&$\\beta_{10}$&$\\beta_{1c}$&$\\beta_{12}$&$\\beta_{20}$&$\\beta_{2c}$&$\\ln(\\tau_1)$&$\\ln(\\tau_2)$ \\\\ ";
        
        print "\\hline";
        print "\\multicolumn{1}{l}{} &  \\multicolumn{7}{c}{\\emph{N=" ""$+ftocv(obshead,0,0) "}} \\\\ \\hline";
        resmat;
        print "\\hline";
        print "\\multicolumn{8}{l}{Kept simulations: " ""$+ftocv(numkeep,0,0)$+"}\\\\";
        output off;    
    elseif j == rows(obsformc);
        output file = "c:/gauss50/research/throbit/mc3.out" on;
        print "\\hline";
        print "\\multicolumn{1}{l}{} &  \\multicolumn{7}{c}{\\emph{N=" ""$+ftocv(obshead,0,0) "}} \\\\ \\hline";
        resmat;
        print "\\hline";
        print "\\multicolumn{8}{l}{Kept simulations: " ""$+ftocv(numkeep,0,0)$+"}\\\\";
        print " \\hline \\hline ";
        print "\\end{tabular}";
        print "\\end{center}";
        print "\\end{table}";
        output off;
    else;
        output file = "c:/gauss50/research/throbit/mc3.out" on;
        
        print "\\hline";
        print "\\multicolumn{1}{l}{} &  \\multicolumn{7}{c}{\\emph{N=" ""$+ftocv(obshead,0,0) "}} \\\\ \\hline";
        resmat;
        print "\\hline";
        print "\\multicolumn{8}{l}{Kept simulations: " ""$+ftocv(numkeep,0,0)$+"}\\\\";
        output off;
    endif;

    /* Print Hard Cases Table */
    
    amper3= "&"$|"&"$|"&";
    slashes3  = "\\\\"$|"\\\\"$|"\\\\";
    caption = "Truth"|"Mean Throbit"|"RMSE";
    if ismiss(resthl)==1;
        meanthl = zeros(7,1);
    else;
        meanthl = meanc(resthl);
    endif;
    b10 = truthth[1]|meanthl[1]|rmseth2[1];
    b1c = truthth[2]|meanthl[2]|rmseth2[2];
    b12 = truthth[3]|meanthl[3]|rmseth2[3];
    b20 = truthth[4]|meanthl[4]|rmseth2[4];
    b2c = truthth[5]|meanthl[5]|rmseth2[5];
    
    lnt1 = truthth[6]|meanthl[6]|rmseth2[6];
    lnt2 = truthth[7]|meanthl[7]|rmseth2[7];
    
    b10s = "" $+ftocv(b10,0,2);
    b1cs = "" $+ftocv(b1c,0,2);
    b12s = "" $+ftocv(b12,0,2);
    b20s = "" $+ftocv(b20,0,2);
    b2cs = "" $+ftocv(b2c,0,2);
    
    lnt1s = "" $+ftocv(lnt1,0,2);
    lnt2s = "" $+ftocv(lnt2,0,2);
    
    resmat2 = caption$~amper3$~b10s$~amper3$~b1cs$~amper3$~b12s$~amper3$~b20s$~amper3$~b2cs$~amper3$~lnt1s$~amper3$~lnt2s$~slashes3;

    if j == 1;    
        output file = "c:/gauss50/research/throbit/hardcase2.out" reset;
        print"\\begin{table}";
        print"\\caption{\\label{hardcase2} Monte Carlo Simulation of Throbit, ``Hard'' Cases}";
        print"\\begin{center}";
        print"\\begin{tabular}{lccccccc}";
        print"\\hline \\hline";
        print"&$\\beta_{10}$&$\\beta_{1c}$&$\\beta_{12}$&$\\beta_{20}$&$\\beta_{2c}$&$\\ln(\\tau_1)$&$\\ln(\\tau_2)$ \\\\ ";
        print "\\hline";
        print "\\multicolumn{1}{l}{} &  \\multicolumn{7}{c}{\\emph{N=" ""$+ftocv(obshead,0,0) "}} \\\\ \\hline";
        resmat2;
        print "\\hline";
        print "\\multicolumn{8}{l}{Number of Hard Cases: " ""$+ftocv(numlose,0,0)$+"}\\\\";
        output off;
    elseif j == rows(obsformc);
        output file = "c:/gauss50/research/throbit/hardcase2.out" on;
        print "\\hline";
        print "\\multicolumn{1}{l}{} &  \\multicolumn{7}{c}{\\emph{N=" ""$+ftocv(obshead,0,0) "}} \\\\ \\hline";
        resmat2;
        print "\\hline";
        print "\\multicolumn{8}{l}{Number of Hard Cases: " ""$+ftocv(numlose,0,0)$+"}\\\\";
        print " \\hline \\hline ";
        print "\\end{tabular}";
        print "\\end{center}";
        print "\\end{table}";
        output off;
    else;
        output file = "c:/gauss50/research/throbit/hardcase2.out" on;
        print "\\hline";
        print "\\multicolumn{1}{l}{} &  \\multicolumn{7}{c}{\\emph{N=" ""$+ftocv(obshead,0,0) "}} \\\\ \\hline";
        resmat2;
        print "\\hline";
        print "\\multicolumn{8}{l}{Number of Hard Cases: " ""$+ftocv(numlose,0,0)$+"}\\\\";
        output off;
    endif;
    if j == 1;
        output file = "c:/gauss50/research/throbit/diag3.out" reset;
    else;
        output file = "c:/gauss50/research/throbit/diag3.out" on;
    endif;
    if j == rows(obsformc);
        print "I have finished with N=" ""$+ftocv(obshead,0,0) ".";
    else;
        print "I have finished with N=" ""$+ftocv(obshead,0,0) ", moving on to " ""$+ftocv(obsformc[j+1],0,0);                
    endif;
    output off;
    j=j+1;
endo;





